Fonte: Storymap UFRRJ
Introdução à Análise Estatística Espacial
O que é Análise Estatística Espacial ?
São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;
Segundo Bailey & Gatrell (1995), “Análise estatística espacial pode ser realizada quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados”
A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?
Origem da Estatística Espacial
- Dr. John Snow (1813-1858) \(\rightarrow\) Considerado pai da Epidemiologia Moderna
Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.
Objetivos da Estatística Espacial
Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)
Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)
Dependência Espacial ou Autocorrelação Espacial
“Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)
“Todas as coisas se parecem, coisas mais próximas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia
Tipologia dos Dados Espaciais
Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.
Noel Cressie (1993) divide a estatı́stica espacial em 3 grandes áreas:
Dados de processos pontuais
Dados de geoestatística
Dados de área
Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados. Exemplo:
Padrões Pontuais
O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.
Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.
Neste caso, o dado aleatório de interesse é a localização espacial do evento.
Estimativas de Kernel (ou Mapas de Calor)
- Estimador de intensidade de distribuição de pontos:
\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004
- Localização da ocorrência de casos de dengue em Belo Horizonte/MG
Geoestatística
Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.
Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas atribuidos a eles uma medida, como por exemplo: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.
A determinação experimental do semivariograma, para cada valor de \(x_i\), considera todos os pares de amostras \(z(x)\) e \(z(x+h)\), separadas pelo vetor distância \(h\), a partir da equação:
\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \]
- Sendo \(\hat{\gamma}(h)\) o semivariograma estimado e \(N(h)\) o número de pares de valores medidos, \(z(x)\) e \(z(x+h)\), separados pelo vetor \(h\).
- Esta fórmula, entretanto, não é robusta. Podem existir situações em que variabilidade local não é constante e se modifica ao longo da área de estudo (heteroscedasticidade).
- Exemplo: Mapa sobre o teor de argila no solo.
Dados de Área
Na análise de áreas o atributo estudado é em geral resultado de uma medida sumáraria (ex: contagem, taxas, médias, etc.) em uma determinada área bem definida, por exemplo um polígono.
Tal polígono cuja forma pode ser complexa bem como as relações de vizinhança.
O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.
Mapa Temático
- Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998
Matriz de Vizinhança
Moran Global, Moran Local e Lisa Map
Moran Global
Moran Local
- Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia
Modelos de Regressão Espacial
Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial
Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);
Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).
Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:
- Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)
\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]
Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.
A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.
- Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)
\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.
Como trabalhar com Análise Estatística Espacial: Algumas ferramentas “0800”
Serão feitas duas aplicações de Estatística Espacial utilização o pacote estatístico R:
Aplicação I: Utilizando a biblioteca tmap para construção de mapas temáticos
library(tmap)
data(World)
tmap_style("classic")
# Desenhando um rápido mapa temático mundial para esperança de vida.
qtm(World, fill = "life_exp")Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr
- Bibliotecas que serão utilizadas:
library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)- Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.
- Primeiro, vamos fazer um gráfico apenas com as geometrias.
- Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:
- Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.
# Entrando com os dados observados na wikipedia
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51, 50, 31, 15,
25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17),
com_rede = c(0.273, 0.412, 0.313, 0.177, 0.513, 0.696, 1.000, 0.974, 0.280, 0.065,
0.191, 0.449, 0.916, 0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353,
0.405, 0.096, 0.400, 0.352, 0.998, 0.347, 0.129))# Construindo o mapa com ggplot
brasil.ufs |>
left_join(acesso_san, by = "code_state") |>
mutate(com_rede = 100*com_rede) |>
ggplot(aes(fill = com_rede), color = "black") +
geom_sf() +
scale_fill_viridis(name = "Municípios com rede de esgoto (%)", direction = -1) +
xlab("Longitude") +
ylab("Latitude") +
annotation_scale(location = "bl") +
annotation_north_arrow(location = "br") +
theme_minimal()- Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).
# Juntando o banco com os dados + malha
coord_pontos <- brasil.ufs |>
left_join(acesso_san, by = "code_state") |>
mutate(com_rede = 100*com_rede) |>
st_centroid()# Construindo o mapa com ggplot
ggplot(brasil.ufs)+
geom_sf() +
geom_sf(data = coord_pontos, aes(size = com_rede), col = "blue", alpha = .65,
show.legend = "point") +
scale_size_continuous(name = "Municípios com rede de esgoto (%)") +
xlab("Longitude") +
ylab("Latitude") +
annotation_scale(location = "bl") +
annotation_north_arrow(location = "br") +
theme_minimal()- Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet
data.frame(st_coordinates(coord_pontos),
com_rede = coord_pontos$com_rede,
UF = coord_pontos$name_state) |>
leaflet() |>
addTiles() |>
addCircleMarkers(~ X, ~ Y,
label = ~ as.character(paste0(UF, ": ", com_rede, "%")),
labelOptions = labelOptions(textsize = "13px"),
radius = ~ com_rede/10,
fillOpacity = 0.5)Aplicação II: Análise exploratória utilizando os dados de dengue em Dourados/MS
Nesta apresentação serão utilizados os dados que fizeram parte da dissertação de Isis Rodrigues Reitman intitulada “SAÚDE E AMBIENTE URBANO: A RELAÇÃO DE INCIDÊNCIA DE DENGUE E AS DISPARIDADES ESPACIAIS EM DOURADOS – MS”, apresentada no Programa de Pós-Graduação em Geografia da Universidade Federal da Grande Douradosos/MS, em abril de 2016.
Essa aplicação também se encontra no Curso de Estudos Ecológicos ministrado para o curso de Pós-Graduação em Epidemiologia em Saúde Pública em 2019, pelos pesquisadores Oswaldo Gonçalves Cruz (PROCC/FIOCRUZ) e Wagner Tassinari (DEMAT/ICE/UFRRJ)
OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)
- Biliotecas do R que serão utilizadas
- Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS
unzip('dados/pop2010.zip', exdir='dados')
pop2010 <- read_csv('dados/pop2010.csv')
unzip('malhas/Setor_UTM_SIRGAS.zip', exdir='malhas')
setor.sf <- read_sf('malhas/Setor_UTM_SIRGAS.shp', crs = 31981)
unzip('malhas/contorno.zip', exdir='malhas')
contorno.sf <- read_sf('malhas/contorno.shp', crs = 31981)OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)
- Fazendo um join com as tabelas com os setores censitários + população
setor.sf <- setor.sf %>% mutate (idsetor = as.numeric(CD_GEOCODI)) %>% inner_join(pop2010,by='idsetor')- Lendo e plotando os casos de dengue georreferenciados em Dourados/MS
unzip('dados/dengue_dourados.zip', exdir='dados')
casos <- read_csv('dados/dengue_dourados.csv')
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)ggplot(setor.sf) +
geom_sf(fill = 'white', color='black') +
geom_sf(data=casos.sf, color='red',size=1) +
theme_void()- Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS
unzip('dados/lixo_dourados.zip', exdir='dados')
lixo <- read_csv2('dados/lixo_dourados.csv')
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)ggplot(setor.sf) +
geom_sf(fill = 'white', color = 'black') +
geom_sf(data=lixo.sf,color='blue',size=1) +
theme_void()Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS
- Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.
ggplot(setor.sf) +
geom_sf(fill = 'white', color = 'black') +
geom_sf(data=lixo2.sf,color='blue',size=1) +
theme_void()- Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)
ggplot(setor.sf) +
geom_sf(aes(fill=pop)) +
geom_sf(data=casos.sf,color='red',size=0.7, aes(colour = "Caso"),
show.legend = "point") +
geom_sf(data=lixo2.sf,color='salmon',size=1, aes(colour = "Lixo"),
show.legend = "point") +
scale_fill_distiller(palette ="PuBu", direction = 1) +
scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
theme_minimal()- Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?
Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.
lixo_buffer <- st_buffer(lixo2.sf, 500)
ggplot(setor.sf) +
geom_sf(aes(fill=pop)) +
geom_sf(data=lixo_buffer,color='gray', fill = "transparent", size=0.4) +
geom_sf(data=casos.sf, color='red',size=0.7) +
scale_fill_distiller(palette ="PuBu", direction = 1) +
scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
theme_minimal()- Represntando os casos e o lixo de forma interativa
tm_shape(setor.sf) + tm_borders("black") +
tm_shape(casos.sf) + tm_dots("red") +
tm_shape(lixo.sf) + tm_dots("green") +
tm_shape(lixo_buffer) + tm_borders("blue") +
tmap_mode("view")- Convertendo o dado de pontos (padrão pontual) para dados de área
## conta casos por setor
casos.sf$contador <- 1
casos <- setor.sf |>
st_join(casos.sf) |>
filter(CLASSI_FIN == 1) %>% ## seleciona somente os casos confirmados
group_by(ID) |>
summarise(casos = sum(contador))
st_geometry(casos) <- NULL ## remove atributos de geometria
## numero de depositos de Lixo por setor
lixo.sf$contador <- 1
lixo <- setor.sf |>
st_join(lixo.sf) |>
group_by(ID) |>
summarise(lixo = sum(contador))
st_geometry(lixo) <- NULL ## remove atributos de geometria
# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |>
left_join(lixo,by = 'ID') |>
left_join(casos,by = 'ID') - Plotando o mapa temático dos casos por setor censitário
- Plotando o mapa temático dos pontos de coleta de lixo por setor censitário
- Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário
setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0 # Transformando os missings em zero
summary(setor.sf$tx) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 1.736 4.065 4.311 56.399
- Plotando a distribuição da incidência em Dourados/MS
library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")
ggplot(setor.sf) +
geom_sf(aes(fill = tx), color = 'black') +
scale_fill_gradientn(colours = pal) +
ggtitle("Taxa de incidência de Dengue") +
theme_void()- Kernel por atributos
Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.
Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA
- Extraindo os centróides dos polígonos em Dourados/MS
centroides <- st_centroid(st_geometry(setor.sf))
# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c('X','Y')
plot(centroides.sp)- Colocando os pontos no formato sp
centroides.ppp <- ppp(centroides.sp$X,centroides.sp$Y, dourados.w)
plot(centroides.ppp,pch=19,cex=0.5)- Fazendo o kernel por atributo da taxa de detecção
kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)- Construindo a matriz de vizinhança para verificar a autocorrelação espacial
Neighbour list object:
Number of regions: 284
Number of nonzero links: 1726
Percentage nonzero weights: 2.139952
Average number of links: 6.077465
- Iremos precisar da coordenadas dos centróides
setor.sp <- as(setor.sf, 'Spatial') # convertendo em formato sp
coord <- coordinates(setor.sp) # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
- Verificando a malha de conectividade da vizinhança de Dourados/MS
viz.sf <- as(nb2lines(viz, coords = coord), 'sf')
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))
# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) +
geom_sf(fill = 'salmon', color = 'white') +
geom_sf(data = viz.sf) +
theme_minimal() +
ggtitle("Vizinhança por \n conectividade") +
ylab("Latitude") +
xlab("Longitude")
mapa.viz- Obtendo a correlação da taxa de incidência de dengue Dourados/MS
Moran I test under randomisation
data: setor.sf$tx
weights: pesos.viz
Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.524720303 -0.003533569 0.001128660
- Plotando o correlograma
Spatial correlogram for setor.sf$tx
method: Moran's I
estimate expectation variance standard deviate Pr(I) two sided
1 (284) 0.52472030 -0.00353357 0.00112866 15.7239 < 2.2e-16
2 (284) 0.23776816 -0.00353357 0.00048540 10.9525 < 2.2e-16
3 (284) 0.09536593 -0.00353357 0.00028358 5.8729 4.282e-09
4 (284) -0.02063378 -0.00353357 0.00019736 -1.2172 0.22351
5 (284) -0.07589158 -0.00353357 0.00016732 -5.5939 2.221e-08
6 (284) -0.06448448 -0.00353357 0.00016165 -4.7940 1.635e-06
7 (284) -0.02708340 -0.00353357 0.00016725 -1.8210 0.06861
8 (284) -0.02744543 -0.00353357 0.00018328 -1.7662 0.07735
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)
5 (284) ***
6 (284) ***
7 (284) .
8 (284) .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.
setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[,5]
tm_shape(setor.sf) +
tm_polygons(col='pval', title="p-valores", breaks=c(0, 0.01, 0.05, 0.10, 1), border.col = "white", palette="-Oranges") +
tm_scale_bar(width = 0.15) - Moran Local (Lisa Map) da taxa de incidência Dourados/MS
Local Morans I Stand. dev. (N) Pr. (N) Saddlepoint Pr. (Sad) Expectation
1 1 -0.009885292 -0.01575255 1.0125682 -0.03600714 0.9712767 -0.003533569
2 2 0.226981101 0.46511580 0.6418485 0.70056077 0.4835772 -0.003533569
3 3 -0.010910944 -0.01829621 1.0145974 -0.04041673 0.9677609 -0.003533569
4 4 0.484675519 1.31014141 0.1901480 1.43930439 0.1500643 -0.003533569
5 5 0.018648578 0.05952726 0.9525322 0.09384037 0.9252360 -0.003533569
Variance Skewness Kurtosis Minimum Maximum omega sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001369880 -0.01569409
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400 0.0028119325 0.44472643
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001590844 -0.01822753
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199 0.0066304485 1.08297547
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199 0.0005595065 0.05929966
sad.u
1 1 -0.01569909
2 2 0.49831660
3 3 -0.01823490
4 4 1.59298211
5 5 0.05942125
[ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[,1]
library(scales)
ggplot(setor.sf) +
geom_sf(aes(fill = MoranLocal), color = 'black') +
scale_fill_gradientn(colours=c("blue", "white", "red"),
values=rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))), guide="colorbar") +
ggtitle("Moran local") +
theme_void()Bibliografia sugerida
Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.
Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995
Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004
Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013
Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.
Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ)
- Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ) link)
Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ)
- Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ) link